clear; clc
%% create model
T1 = 80e-6; T2 = 40e-6;
f0 = 40e3; Rmax = 50e3;

m = Algorithm.BlochEquation(T1, T2, 'solution', 'periodic', 'period', 1/f0);
% m = Algorithm.BlochEquation(T1, T2, 'init', [1 1 0].');
% m = Algorithm.BlochEquation(T1, T2);

%% set magnetic field and pumping 
m.set_Omega('omegaZ', 2*pi*f0) ...
 .set_pumping( m.cosine_pump(Rmax, f0, Rmax) ) ...
 .solve().plot_solution();

%%
fval = 0:1e3:80e3;
demod_fun = @(f) m.set_Omega('omegaZ', 2*pi*f).solve().harmonic('y');
demod = arrayfun(demod_fun, fval);

figure;
plot(fval/1e3, [demod.R; demod.X; demod.Y], '.-' )
xlabel('Frequency (kHz)'); ylabel('demod'); legend('R', 'X', 'Y');
grid on; grid minor;